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ABSTRACT 


During  this  funding  period,  a  semi-classical  macro-kinetic  theory  that  describes  the 
dynamic  behavior  of  carriers  in  a  semiconductor  under  the  influence  of  space-time 
varying  fields  has  been  formulated.  The  macro-kinetic  model  is  considerably  easier  to 
implement  numerically  than  Monte  Carlo  methods  or  those  based  on  the  Boltzmann 
Transport  Equation  (BTE).  Moreover,  the  macro-kinetic  model  requires  orders  of 
magnitude  less  computer  time  to  run.  A  Monte  Carlo  method  has  been  developed  for 
obtaining  the  electron  energy  distribution,  transport  parameters,  and  rate  coefficients 
in  multi-valley  semiconductors.  The  procedure  requires  an  order  of  magnitude  less 
time  than  conventional  Monte  Carlo  techniques. 


I.  OVERVIEW  OF  RESEARCH  PROGRAM 

I.  lNTROlM  CTION 

Midi  IWr  Semiconductor  Switches  (UPSS)  arc 
rcccixi.ie  con.-iJerablc  attention  lor  the  promise  they  oiler 
!o  pulse  p. over  technologv?  "  This  promise  is  founded  on 
the  fact  that  with  semiconducting  materials  as  the 
switching  medium,  it  is  possible  to  obtain:  (a)  high  current 
densities  over  large  areas,  (b)  high  carrier  production  and 
relaxation  rates,  (e)  conductivity  modulation  over  several 
orders  of  magnitude,  and  (d)  high  dielectric  field 
strengths.  UPSS  may  be  developed  that  have  lower 
inductance  and  forward  drop,  higher  rep-ratc,  longer 
lifetime,  and  are  more  compact  than  comparable  gas 
switches.  In  addition,  simple  and  compact  (re-usable) 
opening  switches  and  other  pulse  power  devices  (such  as 
frozen  wave  generators)  seem  to  be  realizable  with 
semiconductor  switching  technology?  Optically  or 
electrically  triggered  bulk  and  junction  devices  and 
externally  controlled  bulk  devices  are  being  investigated  in 
a  variety  of  geometries  and  covering  a  wide  range  of 
parameters.  °  Semiconducting  materials  that  have  been 
considered  arc  Si,  III- Vs,  and  diamond. 

With  AFOSR  support,  our  research  activities  have 
focused  on  the  basic  physics  of  HPSS  devices.  Our  aim 
has  been  to  develop  a  quantitative  understanding  of  the 
role  of  the  various  microscopic  processes,  material 
parameters,  trap  dynamics,  and  space-charge  in  shaping 
the  behavior  of  HPSS.  This  knowledge  is  necessary  for 
guiding  the  scaling  of  the  present  low  power  technology  to 
the  regime  of  interest  in  pulse  power  applications. 
Moreover,  with  sufficient  understanding,  it  may  be 
possible  to  tailor  the  electrical  properties  of 
semiconductor  materials  (beyond  that  of  density 
modulation)  for  specific  applications.  During  this  contract 
period,  we  have  developed  models  for  describing  the 
dynamics  of  the  carriers  under  the  influence  of  space-time 
varying  fields.  In  the  next  section,  a  summary  of  the 
accomplishments  made  during  the  period  of  this  contract 
is  given.  Sections  III  and  IV  are  devoted  to  a  detailed 
discussion  of  the  results  obtained. 

II.  REVIEW  OF  ACHIEVEMENTS  AND 
ACTIVITIES  DURING  THE  CONTRACT 
PERIOD 

During  this  past  funding  period  we  have  accomplished 
the  following: 

(1)  A  semi-classical  macro-kinetic  theory  that 
describes  the  dynamic  behavior  of  carriers  in  a 
semiconductor  under  the  influence  of  space-time  varying 
fields  has  been  formulated.  It  is  essentially  a  "properly 
closed"  set  of  moment  equations.  That  is,  a  macro-kinetic 
distribution  is  introduced  and  the  equation  of  evolution 
for  this  distribution  is  used  to  close  the  moment  equations. 

In  this  fashion,  the  transport  parameters  and  rates  that 
appear  in  the  moment  equations  can  be  determined  from 
first  principles  (in  particular,  without  phenomenological 
assumptions  regarding  the  form  of  distribution).  This  is 


particularly  important  for  describing  high  field  transport 
in  "multi-valley"  semiconductor  materials  such  as  gallium 
arsenide  (GaAs)  (see  next  section). 

A  single  valley  macro-kinetic  model  has  been 
developed  and  compared  to  exact  Monte  Carlo 
simulations  of  carrier  dynamics  in  GaAs  and  the  results 
have  been  found  to  be  in  reasonably  good  agreement. 
The  macro-kinetic  model  is  considerably  easier  to 
implement  numerically  than  Monte  Carlo  methods  or 
those  based  on  the  Boltzmann  Transport  Equation  (BTE). 
Moreover,  the  macro-kinetic  model  requires  orders  of 
magnitude  less  computer  time  to  run.  This  theory  has 
been  published  in  a  paper  that  appeared  in  the  Journal  of 
Applied  Physics  (August  1988)  and  is  presented  in  Section 
III.  This  theory  is  now  being  used  to  develop  a  three 
valley  moment  model  of  carrier  transport  in  multi-valley 
semiconductors. 

A  Monte  Carlo  method  has  been  developed  for 
obtaining  the  electron  energy  distribution,  transport 
parameters,  and  rate  coefficients  in  multi-valley 
semiconductors.  The  procedure  requires  an  order  of 
magnitude  Jess  time  than  conventional  Monte  Carlo 
techniques.  The  technique  has  been  discussed  in  a  paper 
that  appeared  in  the  Journal  of  Applied  Physics  (April 
1988)  and  is  presented  in  Section  IV.  At  present,  a 
number  of  papers  relating  other  results  that  have  been 
obtained  are  in  preparation. 

Published  Papers: 

1.  "Nonequilibrium  Macroscopic  Models  of  Carrier 
Dynamics  in  a  Semiconductor,"  J.  Appl.  Phys.  64, 
1220  (1988)  (with  M.  Cheng  and  C.  Wu). 

2.  "Electron  Energy  Distributions,  Transport 
Parameters,  and  Rate  Coefficients  in  GaAs,"  J. 
Appl.  Phys.  63,  2322  (1988)  (with  M.  Cheng). 

A  number  of  presentations  were  also  made  during  this 
period  at  various  workshops  on  HPSS  sponsored  by  both 
ONR  and  AFOSR.  Abstracts  for  papers  presented  at  the 
spring  meeting  of  the  American  Physical  Society  are  given 
in  Section  V.  Also  during  this  period,  one  student  (M. 
Cheng)  received  his  MS  degree  and  is  currently  working 
on  his  Ph.D.  We  were  invited  to  visit  a  number  of 
government  laboratories  that  have  shown  interest  in  our 
results  (Harry  Diamond  and  EDTL/Fort  Monmouth 
Laboratories). 

1.  M.  Kristiansen  and  W.  Portnoy,  Workshop  on  Solid 
State  Switches  for  Pulsed  Power,  Tamaron, 
Colorado,  January  1983. 

2.  B.  Scnitzky,  Workshop  on  New  Direction  in  Solid 
State  Power  Switches,  Polytechnic  University, 
Farmingdalc,  NY,  1985. 

3.  K.  H.  Schoenbach  and  M.  Weiner,  Workshop  on 
Optically  and  Electron-Beam  Controlled 
Semiconductor  Switches,  Norfolk,  VA,  1988. 
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111.  NONEQUILIBRIUM  MACROSCOPIC 
MODELS  OF  CARRIER  DYNAMICS  IN  A 
SEMICONDUCTOR 
I.  INTRODUCTION 

The  dvnamic  be!. a',  tor  of  free  i  .trnas  m  a  semiconduc¬ 
tor  iiikIlt  i he  influence  ot  a  field  may  he  described  m  t lie 
semi-classical  regime  by  the  tune-dependent  distribution 
function.  /(  K.r ,/ )  t  "  here  k  is  the  cat  ner  momentum,  r  is  its 
position,  and  t  is  time)  !  ’  Given  the  initial  state  of  the  carri¬ 
ers.  the  distribution  function  at  any  other  time  may  he  ob¬ 
tained  from  eithet  the  Boltzmann  transport  equation 
I. BTE ) 5  or  from  Monte  Carlo  simulations.'  When  the 
fields  are  changing  in  space  time  it  is.  in  general,  very  diffi¬ 
cult  to  obtain  the  solution  to  the  BTE.  Moreover,  the  Monte 
Carlo  approach,  although  simple  to  implement,  is  very  tune 
consuming  and  in  some  cases  ( depending  on  the  number  of 
test  particles  used)  prohibitive  Once  the  distribution  func¬ 
tion  is  found,  desired  macroscopic  properties  (which  can  be 
measured  )  can  be  calculated  by  averaging  the  corresponding 
microscopic  properties  over  the  distribution. 1 

An  alternate  approach  for  obtaining  the  macroscopic 
state  of  the  earners  is  in  terms  of  moments  of  the  distribu¬ 
tion.  In  general,  an  exaet  description  of  the  state  requires  an 
infinite  set  of  moments  (this  is  equivalent  to  the  fact  that  we 
need  an  infinite  set  of  moments  to  specify  the  distribution  /). 
These  moments  obey  a  hierarchy  (infinite  set)  of  equations 
obtained  by  taking  moments  of  the  BTE.'  It  will  be  assumed 
that  for  the  situations  of  interest  a  finite  set  of  moment  equa¬ 
tions  can  be  used  to  describe  the  behav  ior  of  the  carriers.'1  To 
assess  the  accuracy  of  this  finite  set,  the  results  for  some 
representative  cases  must  be  compared  with  those  obtained 
from  the  exact  distribution  function. 

This  paper  focuses  on  how  to  obtain  a  closed  set  of  mo¬ 
ment  equations  which  are  valid  in  the  presence  of  space-time 
varying  fields.  This  subject  has  recently  received  consider¬ 
able  attention  because  of  the  deficiency  of  the  drift -diffusion 
equation  in  the  analysis  of  high-frequency  and  submicron 
devices. 10-1  ’  The  drift-diffusion  equation  is  based  on  the  as¬ 
sumption  that  the  carrier  momentum  distribution  is  m  equi¬ 
librium  with  the  local,  instantaneous,  applied  field.  This  re¬ 
stricts  the  validity  of  the  equation  to  slowly  varying  fields 
and  large-size  devices.1'  However,  for  fscnticlassical)  de¬ 
vices  with  very  small  spatial  and/or  temporal  scales,  non¬ 
equilibrium  phenomena,  such  as  velocity  overshoot."  " 
dominate  the  transport  behavior  of  carriers.  In  this  case,  the 
use  of  BTE  or  an  equivalent  nonequilibrium  set  of  moment 
equations  is  essential. 

A  number  of  approaches  have  been  proposed  for  obtain¬ 
ing  .1  closed  set  of  moment  equations.  These  approaches  can 
be  divided  into  three  general  categories,  each  based  on  the 
respective  assumption  that  (a)  for  small  electric  fields  the 
distribution  function  can  be  represented  by  two  terms  in  an 
expansion  in  terms  of  Legendre  polynomials,1"' 16  (b)  the 
distribution  function  is  a  displaced  Maxwellian17-19  when 
the  carrier  concentration  is  sufficiently  high,  and  (c)  the 
unknown  variables  and  coefficients  appearing  in  the  mo¬ 
ment  equations  are  assumed  to  be  functions  of  mean  energy 
only  and  are  obtained  from  phenqmenological  equa¬ 


tions  1  1  The  first  approach  is  closest  in  spirit  to  a  deriva¬ 
tion  bused  on  first  principles.  This  is  important  since  it  can 
prov  tde  guidance  as  to  where  the  model  is  likely  to  fail 

In  the  next  section,  a  mathematical  formulation  of  the 
problem  and  the  foundation  for  an  approach  to  its  solution  is 
presented.  This  approach  is  based  on  the  fact  that  the  mo¬ 
ment  equations  vary  in  a  slower  space/time  scale  than  the 
BTE.  Since  little  is  known  about  the  properties  of  the  scatter¬ 
ing  operators  in  the  BTE  for  the  case  of  a  semiconductor,  we 
have  used  the  physical  knowledge  derived  from  the  moment 
equations  to  implement  an  averaging  scheme  over  the  fast 
variations  of  the  BTE.  It  is  shown  that  it  is  sufficient  to  use 
the  characteristic  times  of  the  moment  equations  to  effect  a 
truncation  scheme  and  prescribe  a  procedure  for  arriving  at 
a  closed  set  of  equations.  This  approach  has  also  been  dis¬ 
cussed  in  connection  with  electron  dynamics  in  gases. ::  In 
Sec.  III.  closed  sets  of  moment  equations  are  derived  for 
three  levels  of  description.  These  equations  have  been  used  to 
describe  the  behavior  of  electrons  in  GaAs  subjected  to  step 
fields.  Although,  the  method  presented  in  this  paper  can  be 
used  to  derive  a  multivalley  macroscopic  model,  a  single¬ 
valley  model  has  been  used  for  GaAs.  The  results  for  the 
average  velocity  and  mean  eneigy  obtained  with  this  formu¬ 
lation  are  compared  with  those  obtained  using  Monte  Carlo 
methods.  For  the  Monte  Carlo  calculations,  a  three-valley 
model  for  GaAs  has  been  used.  The  results  from  these  two 
models  are  in  reasonably  good  agreement.  However,  the  sin¬ 
gle-valley  model  in  this  case  does  not  provide  an  accurate 
description  of  the  behavior  at  electric  fields  above  the  thresh¬ 
old  field  (about  3.5  kV/cm  for  GaAs)  for  the  Gunn  effect 
due  to  strong  intcrvalley  scatterings.  Some  concluding  re¬ 
marks  are  given  in  Sec.  V. 

II.  FORMULATION  OF  THE  NONEQUILIBRIUM 
MACROSCOPIC  DESCRIPTION 

For  simplicity  in  notation  (so  that  subscripts  referring 
to  different  types  of  carriers  need  not  he  introduced),  the 
discussion  will  focus  on  the  dynamics  of  electrons  in  a  single 
valley.  A  similar  treatment  holds  for  electrons  in  other  val¬ 
leys  (also  for  holes  in  the  valence  band ),  with  proper  consid¬ 
erations  given  to  interactions  between  the  various  valleys 
and  types  of  carriers.  Extension  to  multivalley  description  is 
straightforward.  The  situations  of  interest  may  he  described 
by  a  distribution  function  in  (K,r)  space,  /( K.r,/).  This  func¬ 
tion  obeys  the  BT  E;  namely,1 : 

<?,/  +  vVr/+  (<7/yi)E-V„/=  1(f).  (1) 

where  vis  the  macroscopic  carrier  velocity,  K  /f(r.r)  is  the 
electric  field  (either  externally  applied  or  arising  from  space 
charge),  and  I  (f)  is  the  linear  scattering  operator,  v  is  de¬ 
fined  as  the  K-space  gradient  of  the  microscopic  energy  e; 
that  is,  v  =  fi~  'VKe (k).  No  specific  form  for  the  operator  / 
need  be  assumed  at  this  time.  This  operator  describes  a  num¬ 
ber  of  physical  processes  (interactions  between  carriers  and 
lattice)  which  occur  in  different  space-time  scales.  The 
scales  of  interest  are  the  fine-grained  (kinetic),  where 


-2- 


oh. lilacs  occur  in  .1  iv  i i!.::  v.iumn^  nine  iIm.iiicc.  .nul 
the  co.ii  sc-u’i  .lined  .  h  \  Ji  >  >,!\  n.nnic  ,  w  Iici  c  i  h.ilij’cs  occm 
ui  a  m.icraviij'iv  ■>.  ,,!c 

A I  I  lie  m.k  i'onc.'iiu  lc\  A.  I  lie  slate  ol  I  lie  cai  i  ici  s  in  a 
semiconductor  is  gi\  cn  In  I  lie  slate  \  ecloi  ; . >1  fmilc  ihmcii 

sioil )  S  ,  w  hose  coinponeiils  are  inoinenls  of  the  disti  ibii 
lion  That  is.  SN  [in  i  r,/ 1 .  ;  I  .  \  |.wluicw  ioiic 
sponds  It'  a  momenl  of  l lie  disinhiilion  ulncli  mat  Iv  a 
scalar,  vector,  or  tensor. !  and  V  correspi  mis  to  the  highest 
momenl  kept  tn  1 1 1- •  description  I  lie  equations  lot  the  mo¬ 
ments  are  obtained  In  taking  appropriately  weighted  inte¬ 
grals  (in  k  space)  of  Fq.  i  1  The  equations  tor  the  first 
three  moments  [namely,  density  >i{  r j).  mean  energy  f(  r ./). 
and  average  momentum  k(  r  J)  ].  are 

d,n  +  V- (nu)  =  vn.  (2a) 

d,  (ne)  4-  V- (ev)  —  </K  /iu  \ne.  (2b) 

d,  (nic)  +  V-  (kv)  —  <?E//t  =  —  v,„hk,  (2c) 

where  the  bracket  implies  an  average  over  the  distribution,  u 
is  the  average  velocity  [  =  f  v/  cIk  —  fi  '  j  VK  £(  k  )/Jk  | , 
and  v,  vf,  and  v,„  are  the  (space-time  dependent )  effect  ivc- 
carrier-gain.  energy  exchange,  and  momentum  exchange 
frequencies.  These  frequencies  arc  defined  by 


vn  =  J~  /(  /  Wk, 

( 3a ) 

v,  ne  —  j  f(K)/(  /  )</k. 

(  'hi 

—  v,„  UK  —  |  K  I(  J  Jc/k. 

i  3c ) 

Since  it  is  difficult  to  ascribe  physical  significance  to  higher- 
order  moments,  their  equations  of  evolution  are  seldom  writ¬ 
ten. 

Unfortunately,  any  finite  set  of  moment  equations  is  not 
determinate."  For  example,  the  set  of  Eqs.  (2)  contains  un¬ 
known  averages  over  the  distributions  (quantities  in  brack¬ 
ets)  and  unknown  rates  [Eqs.  (3) ).  To  calculate  these  un¬ 
knowns  and  thus  arrive  at  a  determinate  set  of  equations  for 
Ss ,  /  needs  to  be  found.  A  similar  problem  arises  in  classical 
gas  kinetics.2 1:4  and  in  electron  kinetics  in  ionized  gases.22 
In  contrast  to  classical  gas  dynamics  and  in  similarity  with 
ionized  gases,  very  little  is  known  about  the  properties  of 
either  /(  / ),  or  the  operator  {q/fi)  E- V.  —  /( /)  in  Eq.  ( 1 ). 
Because  of  this,  a  more  physical  approach  is  proposed  for 
closing  the  moment  equations.  The  key  to  this  approach  is 
the  use  of  information  from  the  macroscopic  equations  to 
effect  the  truncation.  This  is  outlined  below. 

First,  the  moment  equations  are  ordered  according  to 
their  characteristic  scales.  This  step  requires  a  priori  assump¬ 
tions  about  the  relative  magnitude  of  these  scales.  They  can 
be  made  from  physical  considerations.  In  any  event,  the  or¬ 
dering  that  is  used  needs  to  be  confirmed  after  the  solution 
has  been  found.  Equations  (2a)-(2c)  have  been  ordered  ac¬ 
cording  to  their  characteristic  times.  These  times  are  (in  de¬ 
creasing  magnitude):  r  (effective  carrier  production/loss 
tunc),  r(  (energy  exchange  time),  and  r„,  (momentum  ex¬ 
change  time).  The  higher  moment  equations  would  also 
have  to  be  ordered  accordingly.  It  is  assumed  that  their  char- 


acteiistie  limes  are  smaller  than  those  defined  above.  Note 
that,  m  general.  r(  >  r,„  in  a  semiconductor.11 

Next,  the  number  of  moments  in  the  state  vector  S,v  is 
dciet  muted  I  mm  physical  consideration,  and  from  the  scale 
of  die  desired  description.  Alternatively,  the  number  of  mo¬ 
ments  dial  are  used  determines  the  coarseness  of  the  macro* 
st  ,,j,u  tkset  ipiion.  This  is  because  the  model  is  only  valid  for 
time  scales  of  the  order  of  the  smallest  characteristic  time 
contained  in  a  finite  set  of  equations. 

Finally,  note  that  the  distribution  function  /  which 
satisfies  Eq.  ( 1 )  contains  information  to  all  orders  of  time 
greater  than  a  microscopic  scattering  time.1, 2,25  The  micro¬ 
scopic  scattering  time  is,  in  general,  much  smaller  than  the 
characteristic  times  in  any  finite  set  of  moment  equations. 
The  distribution  function  f  contains  “too  much  informa¬ 
tion"  compared  to  the  state  vector  in  the  scale  of  the 
desired  description.  That  is,  the  unknown  variables  and  pa¬ 
rameters  in  the  SA.  description  do  not  follow  the  fast  varia¬ 
tions  in  /  Thus,  to  obtain  a  determinate  set  of  moment  equa¬ 
tions,  it  is  sufficient  to  use  an  “/’’which  only  contains  infor¬ 
mation  in  the  scale  of  the  moment  equations.  This 
distribution,  the  macroscopic  distribution  function  fM, 
obeys  a  "macroscopic-kinetic  equation.”  The  equation  of 
evolution  for/M,  together  with  the  finite  set  of  moment  equa¬ 
tions,  form  a  closed  set.  This  set  can  be  used  to  describe  the 
nonequilibrium  dynamics  of  the  carriers  in  a  time  scale  cor¬ 
responding  to  the  characteristic  times  of  the  moment  equa¬ 
tions  This  description  is  termed  nonequilibrium  because  fM 
may  be  space-time  dependent.  In  fact,  in  the  time  scale  of  the 
moments,  it  is  equivalent  to  / 

A  number  of  procedures  can  be  used  to  arrive  at  an 
equation  for  fM .  The  objective  in  any  of  these  procedures  is  to 
change  the  scale  of  Eq.  ( 1 )  from  the  microscopic  to  that  of 
the  finite  set  of  moment  equations.26,27  In  this  paper,  the 
technique  proposed  by  Bogoliubov  is  used.26  A  problem 
arises  when  trying  to  solve  the  equation  for /w .  This  has  to  do 
with  the  issue  of  assignment  of  initial  values  to  fM.n  In  this 
paper,  it  will  be  assumed  that  the  moments  of fM  correspond 
to  the  (approximate)  macroscopic  state.  To  solve  the  equa¬ 
tions  for  fM,  it  is  still  necessary  to  know  the  various  micro¬ 
scopic  scattering  processes  for  electrons  in  a  given  semicon¬ 
ductor.  After  solving  for/M,  Eqs.  (2)  and  (3)  can  be  used  to 
describe  the  electron  dynamics  in  the  semiconductor  (see 
Sec.  IV).  The  procedure  outlined  above  is  used  in  the  next 
section  to  obtain  closed  sets  of  moment  equations  valid  in 
three  different  regimes  (time  scales). 


III.  THE  NONEQUIUBRIUM  MACROSCOPIC 
EQUATIONS 

The  approach  outlined  in  the  previous  section  will  now 
be  used  to  obtain  the  nonequiJibrium  macroscopic  equa¬ 
tions.  The  characteristic  times  of  the  macroscopic  equations 
[Eqs.  (2a)-(2c)]  can  be  used  to  define  various  levels  of 
descriptions.  The  most  coarse-grained  description  is  valid 
for  times  in  the  order  of  r  (see  Sec.  II).  From  Eqs.  (2a)- 
(2c),  since  v<vf  <  vm,  there  is  a  time  for  which  the  mean 
energy  and  average  momentum  of  the  carriers  have  relaxed 
to  a  state  of  quasiequilibrium  where  their  subsequent  vari- 


-3- 


ation  is  in  the  scale  of  r,  i.e.,  the  scale  of  the  density  varia¬ 
tions.  For  such  times,  the  macroscopic  evolution  of  the  sys¬ 
tem  can  be  described  in  a  single  time  scale  r.  Thus, 
S,  =  [  n  ( r,f )  ] ;  that  is,  the  macroscopic  state  vector  contains 
a  single  moment,  the  density. 

Progressively  less  coarse-grained  levels  of  description 
can  be  defined  by  systematically  using  an  additional  moment 
in  the  state  vector.  This  assumes  that  the  characteristic  times 
in  the  moment  equations  are  not  degenerate  (i.e.,  character¬ 
istic  times  are  not  equal).  In  this  case,  these  are  times  for 
which  the  higher-order  moments  have  relaxed  to  a  state  in 
which  their  scale  of  variation  is  the  same  as  the  moments 
being  used  in  the  characterization  of  the  macroscopic  state. 
If  there  was  a  degeneracy  (for  example,  v(  1  =  v,„  1 ),  then 
the  corresponding  moments  must  be  collectively  taken  as 
components  of  the  state  vector.  For  cases  of  interest  (carri¬ 
ers  in  a  semiconductor),  the  characteristic  times  are  not.  in 
general,  degenerate.8  Thus,  the  next  less  coarse-grained  level 
of  description  is  in  terms  ofS2  =  [«(r,/),?(r,r)  ] .  This  is  val¬ 
id  for  times  in  the  order  of  '.  From  a  practical  point  of 
view,  the  least  coarse-grained  description  of  interest  is  in 
terms  of  S,  =  [«(r,r),e(r,/),R(r,r)  ],  which  is  valid  for  times 
of  the  order  of  v~ 

To  make  the  equations  that  define  the  macroscopic 
state,  S,,  /=  1,2,3,  determinate,  the  macroscopic-kinetic 
distribution, /M,  must  be  found.  The  time  scale  of/,,  must  be 
consistent  with  the  level  of  description.  Thus,  a  macroscopic 
equation  of  evolution  for  fM  needs  to  be  derived  and  solved 
This  is  carried  out  below  for  each  level  of  description  of  in¬ 
terest;  namely,  S, ,  i  =  1,2,3. 

A.  The  S,  state  (defined  for  times  ~v  ’) 

In  this  case,  only  Eq.  (2a)  and  the  equation  for 
J\t  —f  i,  in  the  r  time  scale  are  necessary  to  describe  the 
evolution  of  the  carriers.  These  two  equations  form  a  closed 
set  The  equation  for  /„  is  obtained  by  changing  the  time 
scale  of  the  BTF.  [Eqs.  ( 1 )  ]  from  the  fine-grained  to  a  r 
scale.  T  his  can  be  achieved  using  a  technique  introduced  by 
Bogoliubov Mathematically,  the  change  can  be  accom¬ 
plished  by  the  following  relation: 

/(K.r.O  =/  [,  [K,«(r,f)].  (4) 

That  is.  in  the  r  scale,  the  space-time  dependence  of  the  dis¬ 
tribution  is  not  explicit,  but  implicit  through  a  dependence 
on  the  density.  The  equation  governing  the  changes  in 
can  he  found  using  Eqs.  (1)  and  (2a).  From  Eq.  (4),  the 
changes  in  /  can  be  written  as 

<7,  /-  <1,  f  1,  d,n,  (5a) 

V, Vr/t.  (5b) 

V./=V,/:M.  (5c) 

Subsequently,  the  subscript  r  is  to  be  dropped  from  the  space 
gradient.  The  time  derivative  of  the  density  may  be  eliminat¬ 
ed  from  Eq  (5a)  by  using  Eq.  (2a),  rewritten  in  the  form 

<),  n  J*  vi/k  Vh  +  J"  /(  f[,)dK.  (6) 

After  placing  Eqs  (5)  and  (6)  into  Eq.  ( 1 ).  the  following 
equation  is  obtained  for  the  distribution: 


t*,  /  (  -  |  <1,  f  V  dK-  Vn  +  J  /( /  'M  )Jk) 

t  v  •  V nc)„  f  -r  ( qVJfi)  ■  VK  f  .1,  =  /(  /  m  )  ■  ( 7 ) 

If  the  deviation  from  spatial  uniformity  is  small,  a  parameter 
(v>  can  be  introduced  into  Eq.  (7)  which  is  indicative  of  this 
assumption.  Using  <S  as  a  basis  for  a  perturbation  expansion, 
the  distribution  may  be  expressed  as 


f  ( KM )  =  y  S'f  v,  ( K ) 

which  substituted  into  Eq.  (7)  results  in  the  following  equa¬ 
tions: 

V  -V  ' '  '  '<1.  /  (-  j  dnf„ydK-Vn} 

t  y  <v  ■  'dj  \,  f  i{  j  )dw 

•■J  J 

k-  V  <V  '  '(d„  f  'u  vV'/i) 

f  y<v|f7,/;,  =  ££'/(/ l,  ).  (8) 

“  n  , 

where  r  —  5  V;  V  — 5V'.  From  Eq.  (8),  the  zeroth  order 
equation  is  found  to  be 

t  ,,/f, )  k ■  v.  /  ( km )  r ,  /( / )  -  a„  f  J  /( /  v,  )dK, 

which  has  the  general  solution 

J  [,  (k .ii)  (9) 

where  j  \,  (  k  )  obeys  the  equation 

(0)K'V,/1,.(k)  =  /(  /'„.,)  (K)  j*  /(  f\,  )dK. 

(  10) 

with  the  condition 

J/'M„(K)rfK  =1.  (II) 

Note  that  the  same  symbol  has  been  used  for  f  (  k.h  )  and 
/  ^  (k).  The  context  in  which  they  are  used  determines  the 
argument.  Equation  ( 10)  has  the  form  of  a  steady  state,  ho¬ 
mogeneous  Boltzmann  equation,  and  /  can  be  identified 
as  the  (zeroth  order)  steady-state  distribution  of  a  homo¬ 
geneous  assembly  of  carriers  in  a  field  defined  by  the  local 
value  of  the  field.  This  is  the  distribution  that  exists  at  (r’.f ) 
if  local  equilibrium  with  the  field  is  assumed.  A  number  of 
techniques  are  available  for  solving  this  equation. '  f  \f  can 
also  be  obtained  using  Monte  Carlo  methods."  Noting  from 
Eq.  (9)  that(9,  / 'w  (K,n)  =  f\,  (k),  the  equation  of 0 ( <5 )  is 
found  from  Eq.  (8)  to  be 

fE-v./i,,  +/'M„  j  /(/'„, )dK 

=  /«,,  j /m,,v  dK-Vn 

-  »/*#„  J  Hf'uJdK,  (12) 

where  the  last  term  is  an  approximation  to  the  correspond- 


mg  term  in  Eq  ;  *  > 

1  inean/uig  Eq  1  !  2  .  and  treating  the  right -hand  stile  as 
a  source.  the  solution  to  Eq.  12'  mav  he  foi  mails  written  as 
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Since  if  dK  —  0  (note  that  \  fd  k  =  n  and  if  dK  has 
been  taken  to  be  one),  the  / 's  above  must  satisfy  the 
conditions  j  fs,  <7k  —  0,  /  =  0.1.  This  implies  that 

i  G  da  =  0. 

The  results  obtaitied  above  can  be  summarized  as  fol¬ 
lows.  The  distribution  function  in  the  r  scale  satisfies  Eq. 
(8).  To  first  order  in  iS,  its  solution  is  given  by 

/.if  =/ i/..  (K)/;(r./)  a-  (S/ \r  ■  V (15) 
\v  here 

J  w  ;  k1  J  Kt  (k)  *  dju  (k). 

This  result  is  the  density  gradient  expansion  which  has 
previously  been  a  priori  assumed  for  the  distribution  func¬ 
tion  '  The  use  of  Eq.  (15)  into  Eq.  (6)  yields  a  diffusion- 
type  equation  fo  '  the  densitv 

d  n  =  —  (\d  —  v,/jf)V7i  +  D  -W'n,  (16) 

where  v,;  v(j{  is  the  effective  drift  velocity,  and  D  is  the 
diffusion  tensor  under  uniform  field  condition  at  the  value  of 
the  local  field.  vdx  is  the  contribution  to  the  drift  velocity 
resulting  from  the  fact  that  J  f<  f)dKf=0.  These  quantities 
are  defined  as 

vj  j  v/w.  d ^  5  j  I(f'„"  )dK, 

D  —  —  b  J"  v/  ^  dK. 

These  coefficients  depend  on  the  applied  field  through  the 
fxl  s.  They  can  be  tabulated  as  a  function  of  field  by  numeri¬ 
cally  solving  for  the/w’s  and  using  the  above  equations. 
Once  these  coefficients  have  been  evaluated,  Eq.  (16)  may 
be  used  to  describe  the  evolution  of  the  carriers.  Thus,  in  the 
time  scale,  the  above  results  correspond  to  nonequilibrium 
diffusion  theory,  v,  v, ,  and  can  be  obtained  at  tj  +  ,  using 
Eq.  I  3).  With  this  information,  Eqs.  (2)  are  determinate  at 


i  ,  and  can  then  be  used  to  determine  the  state  S',  at  /,  ,  , . 

A  more  accurate  expression  than  given  above  for  the 
distribution  in  first  order  (  /  )  can  be  obtained  by  explicit- 

lv  i. iking  into  account  the  space-time  variation  of  the  electric 
field.  I  hat  is,  instead  of  Eq.  (4),  let 

/(  k.r.r)  -  f\t  [k,/i(r.r).  £(r,r)]. 

Following  the  steps  subsequent  to  Eq.  (4),  an  expression  for 
f  \,  is  obtained  that  in  first  order  (/  ^  )  is  proportional  to 
the  space  and  time  derivatives  of  the  electric  field.28  With 
this  approach,  no  specific  form  need  be  assumed  a  priori  for 
the  distribution  function.29'9  Moreover,  a  similar  procedure 
can  also  be  used  to  obtain  more  accurate  expressions  for  the 
distribution  in  the  S2  and  S3  states  (see  below).  This  exten¬ 
sion  of  the  theory  presented  is  not  going  to  be  discussed 
further. 

B.  The  Sj  state  (defined  for  times 

In  this  case,  Eqs.  (2a),  (2b),  and  the  equation  for 
fu  =  /  m  >n  r<  time  scale  are  used  to  describe  the  system. 
The  procedure  for  obtaining  the  equation  for  f\,  follows 
along  the  same  lines  as  the  procedure  for  f'M  (Sec.  Ill  A). 
In  this  case,  the  change  in  scale  is  accomplished  by  the  rela¬ 
tion 

/(K.r.r)  =/i,[K,n(r,r),?(r,f)].  (17) 

That  isvin  the  t(  scale,  the  space-time  dependence  of  the 
distribution  is  taken  to  result  from  changes  in  the  density 
and  mean  energy.  Note  that /  2M  changes  in  two  characteris¬ 
tic  time  scales,  r  and  rt.  From  Eq.  (17),  the  changes  in  / 
can  be  written  as 


d,f=d„  f  l,d,n  +  def  l,d,e, 

(18a) 

Vr/=<//2MVrr  + 

(18b) 

v./=v./2„. 

(18c) 

The  time  derivative  of  the  density  and  mean  energy  may  be 
eliminated  from  Eq.  (18a)  by  using  Eqs.  (6)  and  (2b)  re¬ 
written  in  the  form 

d,e=  n~'  ^  -  J  d„f2Me\dK-Vn  -  J  J-e  f\,e  y  dn-Ve 

+  ?E-J yf2M  dK  +  J eI(f2M)d\Aj. 

(19) 

[  The  term  proportional  to  d,  In  n  is  assumed  to  be  small  and 
hence  has  been  neglected  in  Eq.  (19).]  Note  that  in  this  case, 
the  average  velocity  u  is  obtained  from 

rru  =  J  yf  2M  dK. 

(20) 

After  placing  Eqs.  (6),  ( 19),  and  ( 18)  in  Eq. 
lowing  equation  is  obtained  for  the  distribution 

( 1 ).  the  fo! - 
function: 

d„f2u  (-  ja„/2MvJK  Vn+  j~  I(f\,)d k) 

+  d-ff2Mn-'  ^  -  j* d„f2Me  ydK-Vn 

-  j*  d-ffif€  v  g(k'V?  +  •  j~  yf\i  dK 
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t/i  :  */k  )  •  \  V/(  <>„  f  ■  v  V<  .  >,  :  ,, 

-  ! <,!•:  *  v_  /  /i  /  „ !  :i 

1\>  expedite  tiie  algebra  in  obtaining  a  solution  l*>  1  .q.  i  2  1  ' . 
results  obtained  tor  /  i,  [see  Eq  i  15  i  J  will  be  used  as  a 
guide  That  is.  let  /  \,  be  expressed  to  fust  order  m  spatial 
gradients.  bv 

/  ,,  (K.rt.fc)  -  /  \t  (  k\f  )".  -  !  Ik'.tVI 

•-/  W  (K.f  )  V 11  -  j  (  K.f  Vt 


The  objective  for  the  rest  of  this  section  is  to  arrive  at  the 
lowest-order  solution  for  Higher-order  approximations 
and  the  effect  of  the  gradient  terms  are  to  be  considered  m 
the  future.  Thus,  only  the  first  term  in  this  expansion  will  be 
retained.  The  equation  for  f\,  is  found  to  be 


f\t  /(/m.)^k) 


(  ?E' |  v/;w  d K-f  J  el{f\,  Jt/k-J 

-  (q/fi)E  V„ =  /(/;,.).  1 23) 

This  is  a  nonlinear  equation  for  f\,  .  Neglecting  the  first 
term  (assuming  that  the  effective  carrier  gain  integral  is 
small),  and  using  the  following  definitions: 

u  =  u,=ul(f)=  I  v/^A,  (24a) 


r 

e/(  /  y,  )  c/k. 


Equation  (23)  becomes 

(  t/F  u  *  »’,t)  '7,  /;l,  -  <<?/fi>EVK /;w  =/(/;,.  ). 

(25) 

At  this  level  of  approximation  to  /:m(k ,«,?),  Eqs.  (25), 
( (> ) .  and  (  )9)  form  the  closed  set  of  equations  that  describes 
the  evolution  of  the  system.  /:„  is  made  to  satisfy  the  fol¬ 
lowing  normalization  conditions: 


r 

fit  f<« 


)Jk  ■=  t. 


j  J  \,  dtt  ---  1.  (26b) 

The  solution  to  Eqs.  (25)  and  (  26)  can  be  obtained  as 
follows  Performing  a  change  of  variables  from  (?,k)  to 
( m.c.k  ),  vv here 

K  K,  a,  *  K 

i  a,  iv  a  unit  vector  parallel  to  F  ) . 

2  fa,  (/u  . 
uj  t .  u  . 

K  K  . 

and,  using  the  chain  rule, 


it  it  ] 

i)  ■  ,i  —  •  ,). 

<’(  <),  u 

[  to  zeroth  order  in  <)«.<€)]. 

J  hi  t  . 


I  q  ■  2  5  j  l'ec> wiles 
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/■-'  )  <l ,  /  i,  -/(/;,). 


f  or  (5  ( £  £)/£„,  small  (or  /y,  a  slow  varying 

function  ot  <•< ).  the  lowest-order  solution  to  F.q.  ( 27)  satisfies 
the  following  equation: 

/k  -•  /( /).,  )■  (29) 

This  equation  lias  the  form  of  a  steady-state  BTE  with  <?£rq 
as  the  source.  It  is  equivalent  to  Eq.  (10)  that  arises  in  con¬ 
nection  with  the  description  of  the  S,  state.  Thus,  to  lowest 
order  the  macroscopic-kinetic  distribution  function  obeys  a 
steady-state  RTE  in  an  equivalent  field.  This  observation  has 
a  physical  interpretation.  The  actual  field  E  appears  as  a 
source  term  in  the  moment  equations  and  as  such  causes 
changes  in  the  state  S:.  Since  the  equations  describe  the  evo¬ 
lution  ofS;  in  their  characteristic  timescales,  variations  in  E 
are  "filtered"  by  the  equations;  that  is,  as  far  as  changes  in  S, 
are  concerned.  Thus,  it  is  the  "filtered"  field  which  the  state 
really  “sees."  This  "filtered"  field  is  the  equivalent  field  in 
Eq.  (29). 

Equation  (  29)  can  be  solved  numerically  with  £ey  as  a 
parameter.  By  requiring  that 


-  J  e (k) / r/K, 


a  table  for  ?  vs  £vq  can  be  generated.  Moreover,  using  Eqs. 
(3a)  and  (24),  all  unknown  variables/parameters  in  Eqs 
(2a)  and  (2b)  can  be  determined  as  functionals  of  £ot  or 
equivalently  e.  In  this  fashion,  <•  —  v(e).  r,  -  v((e).  and 
u  =  u,(e).  The  system  of  equations  describing  S,  is  closed 

C.  The  Sj  state  (valid  for  times  ~  v”1 ) 

In  this  case,  Eqs.  (2a)-(2c)  and  the  equation  for 
fM  =/Jf  in  the  rm  time  scale  are  used  in  the  description  of 
the  system.  In  the  spirit  of  Eqs.  (4)  and  ( 17),  the  distribu¬ 
tion  function  is  assumed  to  depend  in  space  and  time  as  fol¬ 
lows: 

/(K,r,/)  =/«  | 

Although  obtaining  an  equation  for  is  straightforward, 

its  solution  is  more  difficult  to  find  than  for  cases  S,  and  S:. 
This  stems  from  the  fact  that,  even  in  zeroth  order,  it  is  an 
equation  in  three  variables  (?,k,k). 

In  this  paper,  instead  of  proceeding  to  find  an  equation 
for  /  and  obtaining  its  solution,  an  approximate  expres¬ 
sion  for  / is  presented.  A  more  complete  theory  (in  es¬ 
sence,  a  more  rigorous  derivation  of  the  expression  present- 
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cd '  will  he  l he  subject  vf.t  future  paper 

In  .ill  eases  Mitci  csi ,  the  iit'iicqmlibr mill  maei t'seopie 
equation',  are  solved  using  mmieneal  uiethi'ds  With  those 
methods,  the  held  /fm  each  time  step  is  taken  to  he  constant. 

1  hits,  the  solution  for  in  a  time  step  corresponds  to  the 
approach  to  equilibrium  of  initial  distribution  in  the  pres¬ 
eme  of  a  held  i  i  e  .  an  .n.oal  value  problem  in  each  lime 
slept  F  rom  the  d- ,et. . stints  in  the  previous  sections,  it  is 
known  that  atV  .  ..  sufficient  time,  the  “initial  distiibution" 
w  ill  evolve  i.. to  r  i  if  the  time  step  is  lone  enough )  finis, 
an  approximate  expression  for  J  can  be  obtained  by  as¬ 
suming  that  this  evolution  can  be  modeled  bv  a  relaxation 
process.  Let  the  macroscopic  state  S  .  and  ;  w  be  knovvn  at 
t:  ne  t  .  when  the  held  is  r t  ).  Note  that  the  field  does  not 
change  in  the  interval  and  that  the  distribution 

f  \,  it)  ( i.e.,  the  distribution  at  time  r  )  is  not  the  equilibri¬ 
um  distribution  for  the  field  /fi  r,/  ).  Thus,  m  spirit  of  the 
above  discussion,  the  distribution  at  time/  -  />  t  where// is  a 
continuous  variable)  can  be  written  as 

f  •/'>-  [  f  )  f\, if  ■  pi  ]<• 

•  f  \i 1 !  ■  / 1 1  lor  a  i  t.  t  h)  i 

where  v„  is  the  average  momentum  exchange  frequency  in 
the  interval  i  /  ./  .  ,  );»■.„  p  j  r.„  t  />)///».  That  is./  \(  ap¬ 
proaches/  exponentially,  due  to  the  relaxation  of  the  fast 
component  (  /  \,  —  f  I  resulting  from  momentum  trans¬ 
fer  at  an  average  rate  v,„  .  The  distribution  at  the  end  of  the 
interval  is  found  by  letting  p  =  A/  in  F.q.  (30).  once 
f  '„  ( /  .  ,  )  is  found.  This  procedure  is  to  be  repeated  at  each 
time  interval.  Note  that  to  implement  this  approach  f 
must  also  be  known  at  each  step  This  problem  has  been 
discussed  in  Sec.  Ill  B. 

A  simple  alternative  to  Eq.  ( 30)  is  to  let  /  ’V(  be  equal  to 
J  in  the  calculation  of  the  unknow  n  v  ariables  in  the  macro¬ 
scopic  Eqs.  (2a)-(2c).  This  substitution  leads  to  the  phe¬ 
nomenological  equations  proposed  by  Shur  "  for  space-inde¬ 
pendent  conditions. 

IV.  EXAMPLE:  THE  RESPONSE  OF  A  HOMOGENEOUS 
CONCENTRATION  OF  ELECTRONS  IN  GaAs  TO  A  STEP 
CHANGE  IN  ELECTRIC  FIELD 

In  this  section,  the  response  of  a  homogeneous  concen¬ 
tration  of  electrons  in  GaAs  to  a  step  change  in  electric  field 
is  investigated  using  the  theory  developed  in  Sec.  III.  For  the 
sake  of  simplicity,  it  is  assumed  that  there  is  no  particle  gain; 
that  is,  //(/-,/)  —  const.  Moreover,  a  single-valley  model  has 
been  used  for  the  band  structure.  As  mentioned  in  Sec.  I,  this 
model  has  severe  limitations,  particularly  at  the  values  of 
livid  under  consideration.  However,  the  objective  of  this  ex¬ 
ample  is  to  illustrate  the  application  of  the  theory  presented 
in  the  previous  section  with  a  "model"  calculation.  The  evo¬ 
lution  of  the  electrons  is  discussed  in  the  context  of  the  S, 
tune  scales  This  ts  dictated  by  the  time  scale  of  the  applied 
field  For  this  example,  Eqs.  (2a)-(2c)  reduce  to 

r>: c  v,(e  —  e„)  -f  quE,  (31a) 

'3,  p  =  —  ym  p  +  qE,  (31b) 


wlieie  c, ,  is  the  mean  energy  corresponding  to  the  lattice 
temperature. />  -  Ek.  and  v  ,  /  —  e,  m,  are  obtained  from  Eqs. 
i  'b )  and  ( 3c  )  with  /  =  /  'w. 

These  equations  have  been  solved  numerically  using  fin¬ 
ite  difference  techniques.'1  At  the  yth  time  step,  the  mean 
energy  ft  /)  and  average  momentum/?!  j)  are  obtained  from 
the  discrete  equations,  given  their  values  and  the  rates  at  the 
previous  time  step y-1.  After  substituting  in  Eqs.  (3b)  and 
!  3c)  for  / [Eq.  (30)],  the  rates  at  the  yth  time  step  are 
obtained  from 


[  AM  I  !  Kale  cnetVicients  for  electrons  in  GaAs  computed  in  the  S:  de- 
scnpii'  -n  as  a  function  of  £“ri| .  or  equivalently,  as  functions  of  e. 


(kV  cm) 

e  (eV) 

(?  —  f„)vj(  10"'  eV/s) 

v„  (f )  ( 10'-’  s  ') 

0  1 

0.319 

0.011  6 

2.292 

0.2 

0.033  6 

0.042  16 

2.684 

0  3 

0.034  75 

0.090  6 

2.796 

0  4 

0.035  7 

0.1566 

2  894 

0.5 

0.036  3 

0  239 

2.933 

0  6 

0.037  3 

0.341  3 

2.95 

0  8 

0.038  8 

0.588 

3.046 

I  (1 

0.040  5 

0.896  4 

3.103 

1 1 

0.041  5 

1.062' 

3.13 

1  2 

0.043 

1.254 

3.187 

1.3 

0.044  4 

1  453 

3.216 

1  4 

0  045  8 

1  663 

3.243 

15 

0.047 

1  886 

3.285 

1  h 

0.047  6 

2  112 

3.317 

1  7 

0.048  1 

2.343 

3.339 

1  8 

0.049  2 

2.343 

3  362 

2.0 

0.055  9 

2.892 

3.517 

2  2 

0062 

3.487 

3.534 

2.5 

0.070 

4.285 

3.638 

5.0 

0  201  7 

9.3 

5.714 

too 

0.3127 

12.53 

12.53 

20  0 

0.370  5 

21  86 

20.39 

30  0 

0.408 

31.98 

26.64 

40  0 

0.441  2 

42.6 

31.92 

500 

0  465  5 

49.5 

39.62 

60  0 

0.497  5 

59.22 

44.25 

70  0 

0.522  9 

68.48 

49.04 

K(j)  =  [k(J~  1)  ~  v’,<OI( y)  ] e  '  +  »-;•”(  y), 

(32) 

where  i  =  e,  m ,  and  the  superscripts  correspond  to  rates  in 
the  S3  scales  (see  Sec.  Ill  C).  The  vj°”s  are  obtained  from 
Eqs.  (3b)  and  (3c),  with  /=/«.  frM  is  obtained  using 
Monte  Carlo  methods.32  The  values  for  the  v(,0”s  obtained 
from  the  Monte  Carlo  model32  are  listed  in  1  .ble  I  as  a  func¬ 
tion  of  mean  energy. 

In  these  calculations,  v  has  been  approximated  by 
vm  (7  —  0-  This  is  the  momentum  exchange  frequency  at  the 
beginning  of  the  interval.  Thus,  with  e(y  -  1),  p(  y  -  1). 
and  v,  (y  -  1 )  given,  the  values  of  e(  j),p{  y),  and  v,  ( j)  are 
obtained  by  solving  Eqs.  (31a),  (31b),  and  (32).  This  pro¬ 
cedure  is  repeated  at  each  time  step. 

The  evolution  of  the  mean  energy  and  average  velocity 
of  the  electrons  in  GaAs  subjected  to  a  step  change  in  elect  ric 
field  is  shown  in  Figs.  1  and  2.  Two  cases  are  shown.  These 
correspond  to  two  different  time  dependencies  of  the  electric 
field  (see  Fig.  3).  The  initial  field  is  kept  constant  for  a  time 
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FIG.  1.  Evolution  of  the  electron  mean  energy  as  a  consequence  of  a  step 
change  in  the  applied  field  (see  Fig.  3 )  .  On  the  figures,  the  large  dot.  solid, 
long  dash,  and  small  dot  lines  correspond  to  the  response  obtained  with  the 
S,S;  andS,  descriptions  and  Monte  Carlo  simulation,  respectively,  (a)  and 
(b )  correspond  to  field  changes  shown  in  Figs.  3(a)  and  3(b).  respectively. 


such  that  the  electrons  have  attained  equilibrium  with  the 
field  by  the  time  the  field  begins  to  change.  Also  shown  in 
Figs  1  and  2  arc  the  results  obtained  using:  (a)  Monte  Carlo 
methods,  with  a  three-valley  model12;  (b)  Eqs.  (31a)  and 
(31b)  with  rates  determined  from  the  S:  state  ( i.e.,  by  letting 
/w  =/ ir );  and  (c)  the  S,  state  approximation. 

The  differences  between  the  results  obtained  with  the  S, 
description  and  the  Monte  Carlo  in  simulation  in  Fig.  2(a) 
arise  primarily  from  the  fact  that  the  S3  description  uses  a 
single-valley  representation.  For  such  high  fields  [see  Fig. 
3(a)]  intervalley  scattering  dominates  the  scattering  pro¬ 
cess  of  electrons  in  GaAs.  It  causes  the  slower  approach  of 
the  average  velocity  (relative  to  the  S,  description)  to  the 
equilibrium  state.  For  lower  fields  [see  Fig.  3(b)],  a  very 
small  fraction  of  the  electrons  gain  sufficient  energy  to  popu¬ 
late  the  upper  valleys  through  intervalley  scattering.  Be¬ 
cause  of  this,  the  response  of  the  carrier  distribution  to  the 
change  in  electric  field  is  not  determined  by  intervalley  scat¬ 
tering  [as  it  is  at  higher  fields;  see  Figs.  2(a)  and  2(b)  ];  and 
thus,  it  is  faster  than  at  higher  fields.  At  lower  fields,  the 
single-valley  nonequilibrium  moment  theory  results  for  the 
mean  energy  are  in  very  good  agreement  with  the  Monte 
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FIG.  2.  Evolution  of  the  average  velocity  of  the  electrons  as  a  consequence 
of  a  step  change  in  the  applied  field  (see  Fig.  3).  The  line  symbols  corre¬ 
spond  to  those  used  in  Fig  1.  (a)  and  (b)  correspond  to  field  changes 
shown  in  Figs.  3(a)  and  3(b),  respectively. 


Carlo  results  [see  Fig.  2(b)].  However,  the  results  for  the 
average  momentum  are  observed  to  relax  slower  than  the 
Monte  Carlo  results.  This  is  also  in  part  due  to  the  fluctu¬ 
ations  in  the  Monte  Carlo  results  at  low  fields  for  the  macro¬ 
scopic  rates  obtained  with  the  S:  description.  Presently,  a 
three-valley  S3  model  is  being  implemented.  The  detailed 
multivalley  effect  on  nonequilibrium  dynamics  of  electrons 
will  be  discussed  in  the  future  paper.  The  evolution  of  the 
system  from  the  initial  equilibrium  stale  ( /)  to  the  final  equi¬ 
librium  state  (T)  is  displayed  in  (e,u)  space  in  Fig.  4.  The 
fast  transient  (nonequilibrium)  behavior  obtained  with  the 
S,  approximation  is  clearly  contrasted  with  those  obtained 
from  the  S,  approximation  (which  in  essence  yields  an  evo¬ 
lution  through  a  series  of  equilibrium  states).  As  expected, 
for  fields  changing  in  time  scales  <-  v,  a  description  in 
terms  of  S,  is  not  satisfactory. 

V.  CONCLUDING  REMARKS 

Nonequilibrium  descriptions  of  the  dynamics  of  elec¬ 
trons  in  a  semiconductor  under  the  influence  of  space-time 
varying  fields  have  been  presented.  These  descriptions  are 
valid  in  different  macroscopic  space-time  scales  which  are 
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FIG.  3.  Time  dependence  of  the  applied  field  corresponding  to  (he  results 
shown  in  Figs.  1  and  2-  (a)  high-field  and  (b)  low-field  cases. 


determined  f'om  the  characteristic  scales  of  the  moment 
equations.  The  results  that  have  been  presented  in  this  paper 
correspond  to  the  lowest-order  solutions  of  these  descrip¬ 
tion  In  the  fastest  scale  (S,),  these  lowest-order  results 
have  been  shown  to  be  in  reasonable  agreement  with  those 
obtained  from  a  kinetic  description.  A  number  of  issues  re¬ 
main  to  be  addressed.  Among  these  issues  are  (a)  the 
(more)  quantitative  description  of  the  S,  state,  (b)  the  rela¬ 
tive  importance  of  higher  order  terms  in  the  expansions  of/„ 
and  of  faster  time  scales  ( S4  or  higher ) ,  and  ( c )  the  relation¬ 
ship  between  a  description  in  terms  of  S,  and  a  modal  de¬ 
composition  of  the  distribution  function.  These  issues  will  be 
discussed  in  a  future  publication. 


Mean  Energy  (ev) 


Mean  Energy  (ev) 


FIG  4.  Phase-space  plot  of  the  evolution  of  the  electron  assembly.  The  let¬ 
ters  I  and  F  correspond  to  the  initial  and  final  states.  The  arrows  indicate  the 
direction  of  evolution.  The  solid  line  with  arrow  corresponds  to  the  response 
in  the  S,  description;  the  line  without  arrows  corresponds  to  theS,  descrip¬ 
tion  (a)  and  (b)  correspond  to  field  changes  shown  in  Figs.  3(a)  and  3(b). 
respectively. 
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IV.  ELECTRON  ENERGY  DISTRIBUTIONS, 
TRANSPORT  PARAMETERS,  AND  RATE 
COEFFICIENTS  IN  GaAs 

I.  INTRODUCTION 

l  ''JvT  L-.-n.un  coiid.iionv  ;hv  helntv  tor  of  an  assembly  of 
.  a.::  -ns  01  holes  m  a  semiconductor.  and  influenced  by 
s;' ..  c-nmc  vary ;ne  fields,  is  described  by  the  tnne-depen- 
:wr.  s,m:el.:ss;cal  density  disinbulion./i  k.r.M.  in  contigu- 
r  ::ivl  momentum  tki  space.1  In  multivalley  semi- 
.  on, tact with  a  different  effective  mass  in  each  valley,  a 
description  of  cm  t  ier  beliav  lor  in  terms  of  a  single  distribu- 
tion  fun.  lion  may  not  he  adequate.1  Moreover,  transport 
parameters  atul  rate  coefficients  that  appear  in  a  single-val¬ 
ley  hydrodynamic  model  do  not  reflect  in  such  cases  the 
dy  namics  of  the  carriers.  1  Ins  is  especially  true  for  values  of 
applied  field  where  intervalley  scattering  is  significant.  In 
these  cases,  the  transport  parameters  not  only  depend  on  the 
transport  properties,  but  also  on  the  rate  coefficients.  For 
example,  drift  velocity,  defined  as  the  time  rate  of  change  of 
the  center  of  mass  of  a  group  of  electrons.-'  will  depend  on  the 
averaged  inter-. alley  scattering  rate.  This  is  due  to  the  fact 
that  carriers  mov  mg  w  ith  the  group  are  scattered  into  valley  s 
with  different  masses  (and  consequently  different  dynam¬ 
ics  i  resulting  in  a  change  in  the  center  of  mass  without  ( nec¬ 
essarily  )  any  transport.  Thus,  in  these  situations,  a  multidis¬ 
tribution  description  of  the  behavior  of  the  carriers  is 
desirable  at  both  the  kinetic  [  /,  ( k.r,/ ),  where  a  denotes  a 
particular  valley]  and  hydrodynamic  (in  terms  of  moments 
of  the/,'s)  levels. 

The  distribution  function,  transport  parameters,  and 
rate  coefficients  in  each  valley  can.  in  principle,  be  obtained 
from  Monte  Carlo  simulation'  "  or  from  solution  of  the 
Boltzmann  transport  equation  (BTE)  by  either  iterative  '" 
or  analytical1"  12  techniques.  At  present,  the  Monte  Carlo 
approach  has  a  number  of  advantages  over  the  BTE  ap¬ 
proach:  it  is  relatively  easy  to  implement  a  six-dimension 
(k.r)  space  simulation:  it  can  be  easily  modified  to  accom¬ 
modate  any  number  of  interactions  between  the  carriers  and 
the  background:  and  it  provides  considerable  physical  in¬ 
sight  into  the  behavior  of  the  carriers,  including  fluctuation 
phenomena. 

In  the  Monte  Carlo  approach,  the  accuracy  of  the  re¬ 
sults  depends  (a)  in  transient  situations,  on  the  number  of 
electrons  used  in  the  simulation,  and  (b)  in  steady  state,  on 
the  total  number  of  scatterings.  A  major  drawback  of  the 
Monte  Carlo  approach  is  that  the  simulation  takes  a  consid¬ 
erable  amount  of  computer  time,  even  when  very  few  elec¬ 
trons  are  used.  This  becomes  more  serious  when  simulating, 
for  example,  the  behavior  of  electron  in  a  multivalley  semi¬ 
conductor  subjected  to  low  electric  field.  In  this  case,  a  very 
small  fraction  of  the  electrons  gain  sufficient  energy  to  popu¬ 
late  the  upper  valleys  through  intervalley  scattering.  For  ex¬ 
ample.  in  GaAs,  less  than  2%  of  the  electrons  are  in  the  X 
valleys  w  hen  E  —  20  kV/cm  ( see  Fig.  4).  Thus,  to  obtain  an 
accurate  representation  of  the  behavior  of  the  electrons  in 
the  X  valleys,  it  is  necessary  to  use  in  the  order  of  tens  of 
thousand  electrons  in  the  Monte  Carlo  calculation.  Such  a 
large  number  of  electrons  make  this  approach,  in  some  cases, 
prohibitive. 

To  simplify  the  computational  aspects  of  the  Monte 
Carlo  approach,  Rees9  1 '  introduced  the  concept  of  self-scat¬ 
tering.  However,  when  the  carrier  scattering  rates  are  in¬ 


creasing  (unctions  of  energy .  standard  implementations  of 
this  concept  leads  to  a  very  large  number  of  fictitious  scatter¬ 
ing  events  along  the  carrier  trajectories,  and  thus  results  in  a 
further  increase  in  the  computation  time.  We  have  devel¬ 
oped  a  technique  for  reducing  the  number  of  self-scattering 
events.  Consequently,  tor  a  given  CPU  (central  processing 
unit )  time,  more  test  particles  can  be  used  in  the  simulation. 
This  results  in  a  reduction  in  the  fluctuation  of  the  calculated 
quantities.  This  technique  is  discussed  in  the  next  section  in 
the  context  of  two  generic  time  dependencies  of  the  applied 
field,  namely,  de  and  a  step  change.  In  Sec.  IV,  results  from 
simulations  of  the  steady-state  behav  ior  of  electrons  in  a 
three-valley  model  of  GaAs  are  presented.  Concluding  re¬ 
marks  are  given  in  Sec.  V. 

II.  THE  MONTE  CARLO  TECHNIQUE 

In  the  Monte  Carlo  approach  for  simulating  the  behav¬ 
ior  of  carriers  in  semiconductors  and  influenced  by  space- 
time  varying  fields,  the  initial  distribution  of  carriers 
[  /(k.r.O)  ]  is  specified.  An  initial  number  of  test  carriers  are 
then  selected  that  are  representative  of  this  distribution,  and 
their  evolution  simulated  using  statistical  methods.2'"1 1  In 
this  paper,  it  is  assumed  that  energy  e  of  an  electron  is  related 
to  the  wave-vector  k  through  the  equation4  14 

h2k2/lm=e(  \  +cxe)  ,  (la) 

where  h  is  the  Planck  constant  divided  by  2~,  in  is  the  effec¬ 
tive  mass  of  the  electron  with  zero  energy  in  the  valley,  and  a 
is  the  nonparabolicity  parameter,  m  and  a  depend  on  the 
valley  in  which  the  electron  is  found.  Equation  (la)  repre¬ 
sents  a  nonparabolic  energy  band  with  spherical  constant 
surfaces  and  a  scalar  effective  mass  m.  For  a  nonparabolic 
energy  band  with  ellipsoidal  constant  energy  surface,  the 
e  —  k  relation  is  given  by 


dfAi+AI), 

2  \m,  m,  ) 


f(k)  f  1  -Forefk)] 


where  m,  and  m,  are  the  longitudinal  and  transverse  compo¬ 
nents  of  the  effective  mass  tensor,151''  and  k ,  and  k,  are  the 
longitudinal  and  transverse  components  of  the  wave  v  ector 
of  the  electron.  For  the  values  of  field  of  interest,  Eq.  ( 1 ) 
represents  a  good  approximation  to  more  accurate  represen¬ 
tations  for  the  band  structure.  Moreover,  the  use  of  Eq.  ( 1 ) 
in  the  trajectory  equations  is  consistent  with  the  formulation 
of  scattering  rates  as  functions  of  energy,  which  assumes  an 
energy-momentum  relation  given  by  Eq.  (  1 ). 

Note  that  it  is  sufficient  to  only  discuss  the  case  of  an 
energy  band  with  spherical  constant  energy  surfaces  [Eq. 
( la)  ].  The  resulting  equations  are  also  valid  for  the  ellipsoi¬ 
dal  case  [Eq.  ( lb)  ]  if  m  is  replaced  with  free-electron  mass, 
and  the  wave  vector  and  electric  field  arc  replaced  with  the 
Herring-Vogt  transformed  values.151'’ 

The  flight  of  an  electron  between  scattering  events  is 
calculated  using  the  equation  of  motion  (h  d,  K  =  ^E)  and 
either  Eqs.  (la)  or  (lb).  The  time  t,  between  scattering 
events  is  determined  from  the  equation2  4 


Ri  =  l-exp(^-J  vT  [f(r)]rfr  j  , 


where  K  is  a  uniformly  distributed  random  number  m  the 
interval  [0,1  ],  and  r,  is  the  total  scattering  rate,  which  is  a 
function  of  the  time-dependent  electron  energy. 

The  integral  m  lit).  (2)  cannot,  in  general,  be  evaluated 
analytically.  To  overcome  tins  dilliculty.  a  fictitious  scatter¬ 
ing  event  is  introduced  such  that  the  "new"  total  scattering 
rate  »•’,  would  be  constant.  '  ' '  This  v',  is  taken  to  be  greater 
or  equal  to  the  minimum  constant  that  makes  r  -,(<-)  posi¬ 
tive  for  all  e  in  the  expression  [see  Fig.  1(a)] 

'■/  =  (  f )  -r  vr(€)  ■  (3) 

vsclt-(e)  is  the  scattering  rate  for  the  fictitious  scattering 
mechanism.  This  process  causes  no  changes  in  the  properties 
of  the  electron  along  the  trajectory.  That  is,  the  state  k'  of  an 
electron  after  a  self-scattering  event  is  taken  to  be  equal  to  its 
state  k  before  the  event.  With  vr  in  Eq.  (2)  replaced  by  v', 
(which  is  constant),  the  integral  in  Eq.  (2)  is  evaluated,  and 
the  duration  of  the  free  flight  ts  is  found.  The  procedure  for 
determination  of  the  scattering  mechanism  and  direction 
has  been  described  in  the  other  papers."1 4  From  Eq.  ( 3 ),  note 
that  the  number  of  self-scattering  events  is  always  much 
greater  than  the  number  of  real  scattering  events.  To  reduce 
the  number  of  self-scattering  events  resulting  from  the  use  of 
Eq.  (3),  a  step-shaped  total  scattering  rate  vf(f)  has  been 
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l-'IG  1  Energy  dependence  of  the  total  scattering  rates.  The  solid  lines  rep¬ 
resent  the  total  scattering  rate  vr  (f  )■  The  dashed  lines  represent  (a)  a  con¬ 
stant  total  scattering  rate  v’r,  (b)  a  step-shaped  total  scattering  rate  vf(e), 
and  (c)  a  total  scattering  rate  (f )  given  by  Eq.  (9a).  is  a  maximum 
electron  energy;  e,  is  a  suitable  boundary  between  two  energy  regions  with 
constant  scattering  rates 


proposed.’  vj  ( e )  is  given  by 


where  c ,  is  a  suitable  boundary  between  two  energy  regions 
with  constant  sealtC'  ing  rates,  namely  v,-,  and  vr2  [see  Fig. 
1(b)  [ .  The  constants  v, ,  and  v ri  are  constant  total  scatter¬ 
ing  rates  which  include  self-scattering.  is  then  obtained 
from  Eq.  (3),  as  discussed  in  Ref.  2. 

The  step-shaped  total  scattering  rate,  vf  (e),  outlined 
above  does  not  significantly  decrease  the  simulation  time 
unless  more  steps  are  used.  On  the  other  hand,  the  use  of 
more  steps  cause  difficulty  in  the  implementation  of  the 
scheme.  This  is  because  situations  in  which  an  electron  tra¬ 
vels  across  two  or  more  energy  regions  without  suffering  any 
scattering  have  to  be  considered. 

To  significantly  reduce  the  number  of  self-scattering 
events  while  keeping  the  implementation  of  the  scheme  rela¬ 
tively  simple,  we  propose  the  following  scheme:  (a)  change 
the  integration  variable  in  Eq.  (2)  from  time  (  to  momentum 
k;  (b)  use  a  quadratic  polynomial  to  represent  the  total  scat¬ 
tering  rate  (including  self-scattering)  and  (c)  use  the 
energy-momentum  relation,  Eq.  (1),  to  carry  out  the  inte¬ 
gration. 

Following  this  procedure,  let  vT(e)  in  Eq.  (2)  be  given 
by  Vj-;  that  is, 

v’jTf)  =A  +  Be{\  +ae)  ,  (4a) 

and  using  Eq.  ( 1 ), 

v«r(*)  =A  +B(h2k2/lm)  ,  (4b) 

where  A  and  B  are  constants,  which  are  chosen  from  the 
requirement  that  v5e,r(e )  in  Eq.  (3)  must  be  positive  and  as 
small  as  possible  in  the  energy  range  of  interest  [see  Fig. 
1(c)].  After  changing  the  variables  of  integration  from  t  to 
kz  and  utilizing  Eq.  (4)  and  the  equation  of  motion,  Eq.  (2) 
becomes 


Ah  f  '  —  dkz=  —  In ( 1  —/?,), 

-L  qE  2m  aE 


where  k ^  is  the  initial  value  of  the  longitudinal  component 
of  the  wave  vector  k  ,k2  =  k  \  +  k1  +  k],  and  we  have  tak¬ 
en  E  along  the  z  direction.  E  in  Eq.  ( 5 )  may  be  either  space- 
time  dependent  or  constant.  The  piecewise  application  of 
Eq.  (4b)  in  the  energy  interval  of  interest  results  in  further 
reduction  of  the  number  of  self-scattering  events.  The  result¬ 
ing  total  rate  v"j  is  further  discussed  in  the  next  section. 

Since  space-time  simulations  are  implemented  as  a  se¬ 
ries  of  time  intervals  in  which  the  field  is  constant  in  time,617 
it  is  only  necessary  to  investigate  the  generic  problems  of 
constant  and  step  field  variations.  These  two  cases  are  con¬ 
sidered  below. 

1.  Static  electric  field:  In  this  case,  Eq.  (5)  reduces  to 

‘■+3(*'+! 

~{j^kzo  +  3*'**° +*rf)=o’  (6) 
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where  A.  x  A  ,  -  A  is  the  transverse  component  of  the 
w.oe  wcior  which  is  constant  during  the  free  (light.  Since 
Eq  !  ti  i  hasonlv  one  real  root.  it  simplifies  the  determination 
of  A  at  the  end  ofthe  free  (light  for  a  given  random  number 
K  .  I'lte  duration  of  the  free  flight  /,  is  then  obtained  from  the 
equation  of  motion. 

1  Step  electric  jidd:  The  ease  of  a  step  v  ariation  in  time 
of  the  applied  field  can  be  treated  in  a  similar  way  as  for  a 
constant  held.  Let  the  applied  field  be  in  thee  direction,  and 
equal  to  E,  for  t  i ,  and  E„  for  t  >  tr  To  determine  the  state 
of  the  electron  at  the  end  of  a  free  path,  it  is  necessary  to 
know  t  he  location  of  the  free  flight  relative  to  that  ofthe  step. 
This  gives  rise  to  two  situations: 

t  i )  An  electron  ends  its  free  flight  at  t  t,  or  starts  the 
free  flight  at  t  t,.  In  these  cases  Eq.  (b)  is  solved  for  the 
longitudinal  component  ofthe  wave  vector  A.  at  t  he  end  of 
the  free  path,  within  E  --  E,  or  E  =  Eu .  respectively. 

( li )  An  electron  begins  its  free  flight  at  t  <  /,  and  ends  it 
at  t  >  f,.  In  this  case,  the  integral  in  Eq.  (5)  can  be  decom¬ 
posed  into  two  regions  (  namely,  E  —  E,  and  E  ~  E,, )  and 
integrated  to  giv  e 

k  ~  3(A' ;  +  Tni)k:  +  ^r^(ln( 1  -  *.  >  + 

-  +3*;X.  +  *.u)  =  0.  (7) 

where 

C,  =  {Ah  /qE,  ){k. ,  -  )  =  A(t,  -  t„)  , 

C,  =  ( Bh  72mqE,  )  [  A ;  ( A.,  -  k„, ) 

T  {(A  !|  —  A  )  ]  , 

r„  is  the  initial  time  of  the  free  flight,  and  A.,  is  the  longitudi¬ 
nal  component  of  the  wave  vector  of  the  electron  at  t  = 
The  longitudinal  component  of  the  wave  vector  A.  at  the  end 
of  the  free  flight  is  determined  by  solving  this  cubic  equation. 

111.  ELECTRON  ENERGY  DISTRIBUTIONS,  TRANSPORT 
PARAMETERS,  AND  RATE  COEFFICIENTS  IN  GaAs 

The  Monte  Carlo  technique  presented  in  the  previous 
section  has  been  implemented  to  obtain  the  steady  state  and 
step  response  of  electrons  in  GaAs.  As  mentioned  in  the 
previous  section,  these  are  the  two  generic  problems  which 
form  the  basis  for  simulations  with  arbitrary  time  depen¬ 
dence  of  the  field.  In  this  paper,  the  steady-state  results  are 
discussed.  The  transient  results  are  to  be  discussed  in  a  fu¬ 
ture  publication.  In  Sec.  Ill  A,  further  computational  details 
are  given.  In  Sec.  Ill  B,  the  gain  in  computation  time  that 
has  been  achieved  with  the  technique  presented  in  Sec.  II  is 
illustrated.  The  results  for  the  distribution  function,  trans¬ 
port  parameters,  and  rates  coefficients  are  presented  in  Sec. 
Ill  C. 

A.  Further  computational  details 

In  the  Monte  Carlo  simulations  discussed  in  this  paper, 
a  three-valley  model  ( T,LJ()  of  GaAs  has  been  used.  The 
scattering  processes  that  have  been  taken  into  account  are 
polar  optical,  acoustic  phonon,  intervalley,  and  intravalley 


scatterings. :  ‘  flic  intervalley  separations  Ae,,  and  Ac,  v 
have  been  taken  to  be  equal  to  0.33  and  0.522  eV.  respective¬ 
ly.1'  A  list  ofthe  material  parameters  that  have  been  used  in 
the  simulations  is  given  in  Table  I.  The  values  given  in  paren¬ 
theses  are  those  obtained  by  Pozela  and  Reklaitis1 '  from  best 
fit  to  the  data  of  Houston  and  Evans.'" 

To  get  an  accurate  representation  of  the  steady  state  of 
electron  behavior  in  each  valley,  up  to  80  000  electrons  have 
been  used  in  the  simulations''1  (each  electron  suffers  at  least 
100  collisions).  The  sampling  procedure  has  been  discussed 
in  detail  elsewhere/’ 

B.  Illustration  of  gain  in  computation  time 

To  illustrate  the  gain  in  computational  speed  that  has 
been  achieved  with  the  technique  presented,  we  have  carried 
out  comparative  simulations  with  the  three  total  scattering 
rates  discussed  in  Sec.  II,  namely,  v', ,  v} ,  r") ,  and  v^'.  The 
rate  r'J'  is  a  modification  of  v/  obtained  by  piecewise  applica¬ 
tion  of  Eq.  (4b).  We  have  used  for  the  rate  in  the  T  valley 
in  one  of  our  illustrations.  Note  ’hat  approaches  v)  in  a 
given  interval  when  B  —  0.  Thus,  v’/  is  the  most  general  rate 
that  can  be  used  to  represent  the  total  scattering  rate  while 
still  being  able  to  carry  out  the  integration  in  Eq.  (2).  This 
can  be  considered  as  the  use  of  quadratic  splines  for  the  rep¬ 
resentation  of  the  actual  scattering  rate.  That  is,  in  the  T 
valley,  the  energy  interval  is  divided  into  two  domains 
(e  <  SerL  and  e>  SerL  )  and  Eq.  (4b)  is  applied  to  each 
region.  For  e  <Aert,  B  is  taken  to  be  zero  for  illustration 
( see  below  for  further  comments ) .  The  results  are  shown  in 
Table  II.  For  the  same  number  of  electrons  and  simulation 
time  ( 10  ps),  the  number  of  real-scattering  events  are  ap¬ 
proximately  equal.  The  number  of  self-scatterings  (and  total 
scatterings),  however,  is  considerably  different.  This  also 
applies  to  the  computation  (CPU )  time.  As  seen  from  Table 
II,  the  computation  time  for  the  Monte  Carlo  simulations 
using  t'j  (constant  total  scattering)  is  approximately  an  or¬ 
der  of  magnitude  and  four  times  larger  than  the  computation 
times  for  simulations  using  \\  [Eq.  (4)  ]  and  vf  (piecewise 
constant  total  scattering),  respectively.  Further  reductions 
in  the  number  of  self-scatterings  were  obtained  with  (see 
Table  II). 

At  10  kV/cm,  the  mean  energy  of  electrons  in  the  T 
valley  of  GaAs  is  about  0.23  eV  below  the  threshold  energy 
Acrt.  As  mentioned  in  Sec.  II,  for  this  mean  energy,  the 
electron  scattering  rate  in  T  valley  of  GaAs  is  nearly  con¬ 
stant  and  results  in  a  large  self-scattering  rate  with  the  appli¬ 
cation  of  v\  for  the  whole  interval.  Thus,  in  Table  II,  the 
results  for  v^-(e)  show  approximately  200  self-scattering 
events  (more  than  50%  of  total  scatterings).  This  number  is 
reduced  to  less  than  100  self-scattering  events  with  the  use  of 
v?  as  shown  in  Table  II  (approximately  30%  of  total  scatter¬ 
ings) .  The  number  of  self-scattering  events  is  expected  to  be 
less  if  the  constants  A  and  B  used  for  v£(e )  in  each  energy 
region  are  optimized  by  properly  partitioning  the  energy  in¬ 
terval.  By  dividing  the  total  energy  interval  into  more  sec¬ 
tions  and  applying  to  each  section  an  equation  of  the  form  of 
Eq.  (4),  the  coefficient  At  and  Bt  (where  j  denotes  the yth 
section)  can  be  chosen  to  minimize  the  number  of  self-scat¬ 
tering. 
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TABLE  I.  GaAs  material  parameters  ' 


Density  (g/cm1) 

Sound  velocity  (cm  s) 
Static  dielectric  constant 
Optical  dielectric  const  ant 
1.0  phonon  energy  i  eV  I 
Energy  separation  (eV  i 

r-L 

r.v 


s  In  (S  .17) 

5  74  ■  10'  (S  ’  .  10') 
1 2.0 
io  <i; 

lit)  15  .)(,  (O  l)l(,7  } 

0  VI 

0  522  (0.52) 


EtOOO) 


/.(111) 


AS  ( 100) 


Nonparabolicily  (eV  " 1 ) 

Uhl  (0.62) 

0.461  (0.5) 

0.204  (0.3) 

Effective  mass  (m*/m„) 

0  062 

0.222  (0.17) 

0.58 

Acoustic  deformation 

potential  (cV) 

7.0 

9.2  (7.0) 

9.7  (7.0) 

Intervalley  phonon 

energy  (eV) 

r  (Ocx)  > 

/.(Ml) 

AT  100) 

r(000) 

0 

0.0278  (0.0299) 

0.0299 

/.(HI) 

0  0278  (0  0299) 

0.029  (0.0299) 

0.0293  (0.0299) 

X(IOO) 

0.0299 

0.029.1  (0.0299) 

0.0299 

Intervalley  coupling 

constant  ( 10’  eV/cm) 

r(ooo) 

/.(111) 

X(100) 

r(000) 

0 

1(0.18) 

1 

22(111) 

1(0.18) 

1(0.5) 

05(0.1) 

X( 100) 

i 

0.5(0. 1) 

0.7(1) 

‘See  Refs.  18  and  19. 


C.  Steady-state  behavior  of  electrons  in  GaAs 
1.  Energy  distribution  functions 

To  elucidate  the  physics  of  the  steady-state  behavior  of 
electrons  in  a  three-vallev  model  of  GaAs  we  have  obtained 
the  energy  distribution  of  the  electrons  in  each  valley.  These 
distributions  are  shown  in  Fig.  2. 

At  low  and  medium  fields  ( below  20  kV /cm )  where  the 
population  of  electrons  in  the  X  valleys  is  not  significant, 
these  distributions  have  similar  features  as  those  obtained  by 
Fawcett  and  co-workers4  ”  ' '  and  Conwell  and  Vassell. 14  As 
the  field  increases  from  zero,  electrons  are  heated  up  rapidly 
due  to  invariant  polar  optical  phonon  scattering4  and  the 
distribution  begins  to  flatten.  However,  because  of  strong 
intervalley  scattering,  electrons  with  energy  above  the 
threshold  for  scattering  into  the  L  valleys  ( AerL  )  are  driven 
into  equilibrium  with  the  L  distribution.  This  distribution  is 
nearly  a  Maxwellian  at  the  lattice  temperature.  For  fields 


above  10  kV/cm,  a  population  inversion  is  observed  in  the  T 
valley  in  agreement  with  the  results  of  Fawcett  and  co- 
workers4”-24  and  Conwell  and  Vassell.15  This  inversion  is 
due  to  the  fact  that  intervalley  scattering  is  nearly  isotropic, 
whereas  polar  optical  scattering  is  primarily  forward.  On 
average,  half  of  the  electrons  scattered  from  the  L  to  the  F 
valley  near  the  threshold  energy  Aert  lose  energy  to  the  field 
and  thus  move  to  energy  states  below  b.erL .  Since  intervalley 
scattering  in  this  range  is  zero,  these  electrons  represent  an 
“uncompensated”  source  into  these  states,  thus  creating  the 
observed  population  inversion.  This  type  of  distribution  is 
not  stable  and  should  lead  to  collective  excitations  (plas- 
mons).24  These  excitations,  however,  are  not  taken  into  ac¬ 
count  in  this  model.  However,  in  contrast  to  the  two-valley 
model,  the  population  inversion  is  not  very  pronounced  so 
that  the  plasmon  excitation  rate  should  be  small. 

At  higher  fields  (above  20  kV/cm),  electrons  in  the  L 
valleys  (as  well  as  the  tail  of  the  T  distribution)  begin  to  heat 


TABLE  tl.  Comparison  of  computation  time  for  the  techniques  described  in  Sec.  II.  Simulation  time  is  10  ps.  E  =  10  kV/cm. 


Number  of 

Number  of 

Number  of 

CPU 

Function  for 

total  scattering 

real  scattering 

electrons 

time* 

total  scattering  rate 

for  single  electron 

for  single  electron 

in  simulation 

(min) 

Constant  total 
scattering  rate  v'T 

II  000 

150-190 

2000 

196 

Three-level  step-shaped 
total  scattering  vf 

2300 

- 

150-190 

2000 

49 

Eq.  (4)  rate 
vt 

350-390 

150-190 

2000 

21 

v? 

230-290 

150-190 

2000 

17 

•See  Ref.  21. 
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up  and  the  /.  distribution  i  and  tail  of  T  distribution)  begin 
to  tl.it ten  i  see  f  ig.  2).  Again,  strong  mtervalley  scattering 
drives  the  tail  of  the  L.  distribution  (and  ihe  part  of  the  f 
distribution  above  the  threshold  for  scattering  into  the  .V 
v  ailev  s  :  in  to  equilibrium  with  the  .V  distribution.  This  distri¬ 
bution  is  near]}  a  Maxwellian  at  the  lattice  temperature. 
Tins  behavior  is  similar  to  what  happens  between  the  L  and 
I  \  ailev  s  tit  low  fields. 

At  much  higher  fields  i  above  50  kV/cm ).  the  .V  distri¬ 
bution  also  heats  up.  This  leads  to  the  heating  of  which  in 
turn  heats  up  the  tail  of  the  F  and  the  L  distributions,  and  all 
three  distributions  begin  to  flatten  (see  Fig.  2).  In  this  mod¬ 
el.  there  are  no  other  mechanisms  for  cooling  the  tails  of  the 
distributions.  As  previously  mentioned,  scattering  into  the X 
valleys  prevents  the  population  inversion  in  the  T  valley  to 
become  as  pronounced  as  in  a  two-valley  model.  For  the 
fields  investigated  no  population  inversion  is  observed  in  the 
/.  or  T  valleys  due  to  the  flux  of  electrons  from  the.V  valleys 
near  threshold,  although  the  trend  for  such  a  condition  is 
evident. 

2.  Rate  coefficients 

Macroscopic  (moment)  descriptions  of  carrier  dynam¬ 
ics  in  semiconductors  require  knowledge  of  the  rate  coeffi¬ 
cients  that  appear  in  the  corresponding  equations.  In  a  one- 
moment  description  (in  terms  of  the  continuity  equations 
for  carrier  densities),  the  necessary  rates  are  the  averaged 
carnet  gam  or  loss  rates.  At  applied  fields  for  which  the 
population  ol  the  upper  v  alleys  becomes  significant,  a  niulti- 
v  alley  macroscopic  description  is  desirable.  In  this  case,  car¬ 
rier  gain  or  loss  m  each  valley  is  due  ;  at  lo  arrier  densities 
and  at  fields  for  which  impact  ionization  is  negligible!  to 
inters  alley  scattering. 

With  the  code  described  in  Sec.  II.  we  have  obtained  the 
averaged  (over  the  distribution)  intervalley  scattering  rates 
for  carrier  gain  or  loss  as  a  function  of  applied  field.  The  rates 
are  shown  in  Fig.  3.  These  rates  are  defined  as  follows: 


ELECTRIC  FIELD(kV'cm) 


FIG.  3.  Averaged  transition  rates  as  a  function  of  applied  electric  field.  Sol¬ 
id  lines  represent  the  transition  rates  from  upper  to  lower  valleys.  Dashed 
lines  represent  the  transition  rates  from  lower  to  upper  valleys. 

-  f  S,J(k)/i(k)e/k  , 

where  /  and  j  represent  the  T,  L,  or  X  valleys,  StJ  is  the 
microscopic  intervalley  scattering  rate  from  the  /  valley  to 
the  j  valley,/,'  is  the  distribution  in  the  i  valley,  and  X,  is  the 
population  of  the  i  valley.  Their  field  dependence  follows 
from  the  behavior  of  the  distributions  and  the  microscopic 
intervalley  scattering  rates.  Since  the  microscopic  scattering 
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FIG.  2.  Energy  distributions  of  electrons  in  each 
valley  of  GaAs  for  different  electric  fields.  The 
numbers  correspond  to  the  applied  electric  fields  in 
kV/cm. 
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FRACTION  OF  ELECTRONS  IN  EACH  VALLEY 


rates  fi > 'in  ! !ic  upper  to  lower  valleys  have  no  strong  energy 
dependence. :  a.uj  the  distribution  function  in  the  upper  val¬ 
ley  ^  d«v>  n*'t  eii.uige  sigmtieantly  with  field  (see  Fig.  2),  the 
voiding  of  the  two  fimetions  lead  to  nearly  constant  average 
macroNeopiv  rate1-.  On  the  other  hand,  scattering  from  lower 
to  upper  \  alley  s  wnohcs  the  tail  of  the  corresponding  distri¬ 
butions  since  only  electrons  with  energy  above  the  thresh¬ 
old  can  participate  t .  Since  the  tail  oft  lie  distribution  are  very 
sensitive  to  the  field,  the  folding  of  the  two  functions  is 
strongly  field  dependent. 

At  low  fields  the  average  rates  for  scattering  into  the 
upper  valleys  are  smaller  than  their  respective  inverse  rates. 
As  the  field  increases  above  3.5  kV/cm  (see  Fig.  3),  the 
average  scattering  rate  into  the  L  valleys,  vri ,  becomes  com- 
parable  to  its  inverse  vLr  (higher  than  10%  of  vLr  ),  and  the 
population  of  electrons  in  the  L  valleys  increases  drastically 
(sec  Fig.  4).  Above  35  kV/cm,  the  average  scattering  rates 
into  the  .V  valleys.  r,-.v  and  v,  v,  become  comparable  to  their 
respective  inverses.  vu  and  vv/  .  resulting  in  an  increase  in 
the  population  of  the  .V  valleys  with  a  concomitant  decrease 
in  the  population  of  the  L  valleys  ( see  Fig.  4).  As  previously 
mentioned  at  these  fields,  the  tail  of  the  distributions  heat  up 
due  to  the  heating  up  of  the  population  (see  Fig.  2).  The 
scattering  out  of  the  X  valleys  into  the  L  and  T  valleys  is 
greater  that  out  of  the  L  valleys  into  the  V  valley.  This  has  a 
considerable  effect  on  the  dynamics  of  the  carriers  at  these 
fields  in  steady  state  and  in  transient  situations.  This  is  illus¬ 
trated  in  the  tie**  section  with  regards  to  the  steady-state 
transport  parameters. 

3.  A  verage  velocity  and  mean  energy 

The  technique  discussed  in  Sec.  II  has  permitted  the  use 
of  a  large  number  of  electrons  in  our  simulations.  As  a  conse¬ 
quence,  we  have  been  able  to  observe  steady-state  transport 


properties  which  arc  not  clearly  exhibited  in  simulations 
with  large  fluctuations  i  small  number  of  particles  1.  Here, 
vve  present  the  results  vve  have  obtained  for  the  average  ve¬ 
locity  and  the  mean  energy  for  electrons  in  each  valley. 

The  average  velocity  in  each  valley  is  shown  in  Fig.  5(a) 
as  a  function  of  field.  Also  shown  is  the  total  average  velocity 
<r>  i  averaged  over  all  electrons  irrespective  of  their  valley) 
and  representative  experimental  results.-"-'-'’  The  total 
average  velocity  has  a  maximum  at  a  field  equal  to  3.8  kV/ 
cm  and  rapidly  decreases  with  increasing  field  until  it  be¬ 
comes  nearly  constant  for  fields  above  10  kV/cni  [see  Fig. 
5(a)  j .  This  behavior  is  well  known  and  is  the  result  of  inter- 
valley  scattering  (which  is  isotropic)  and  the  fact  that  the 
electrons  have  a  higher  effective  mass  in  the  upper  valleys. 

At  low  fields  (below  10  kV/cm),  our  calculations  of 
average  velocity  are  in  good  agreement  with  the  results  of 
Ruch  and  Kino.3'  However,  the  calculated  peak  velocity  is 

slightly  lower  than  the  experimental  value  which  occurs  at 
3.5  kV/cm.  In  the  high  field  range  (above  20  kV/cm  ),  our 
calculations  of  average  velocity  fit  the  data  of  Riginos'"  "'un¬ 
well  for  values  of  fields  up  to  70  kV/cm.  These  values  are 
16%  higher  than  those  obtained  by  Houston  and  Evans. 
These  differences  are  reflected  in  the  values  of  the  material 
parameters  used  in  the  calculations  (see  Table  I ).  The  frac¬ 
tion  of  the  electron  population  in  (he  upper  valleys  as  a  func¬ 
tion  of  field  is  shown  in  Fig.  5(b).  This  figure  covers  the 
range  of  fields  near  the  maximum  of  the  total  average  veloc¬ 
ity  and  illustrates  therWrensc  in  (f)  as  the  population  in  the 
upper  valleys  increases.  The  T  valley  average  velocity  <ur  > 
is  observed  to  reach  a  maximum  at  a  larger  field,  namely  5.5 
kV/cm.37  Between  5.5  and  10  kV/cm,  (/»,-  )  makes  a  "transi¬ 
tion”  to  a  lower  value  and  remains  nearly  constant  with  field 
until  at  35  kV/cm  a  second  “transition”  to  a  lower  value  is 
observed.  These  transitions  have,  in  part,  the  same  origin; 
namely,  the  fact  that  intervalley  scattering  is  isotropic.  For 


MG.  4.  Fraction  of  electrons  in  each  valley  as  a 
function  of  applied  electric  field. 
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AVERAGE  VELOCITY  (10  cm/sec) 


fields  above  i lie  first  critical  field,  scattering  from  /.  valleys 
to  F  vallcv  tandomi/es  the  electron  momentum  tit  the  ! 
'■alley  and  tints  decreases  ■ ;  ,  ).  At  high  fields  above  the  sec¬ 
ond  et ttteal  field,  the  population  in  the  /.  valleys  is  observer! 
to  decrease  vv  nil  a  er'ncr'tmtant  increase  tit  the  population  of 
the  A  vallcv  s  ,  see  Fig.  4'.  Since  tiie  scattering  from  the  X 
valleys  into  the  T  valley  is  always  larger  titan  scattering  from 
the  L  valleys,  there  is  an  increase  in  the  number  of  electrons 
with  random  momentum  transferred  into  the  T  valley  above 
45  kY/cm.  This  effect  lowers  the  average  velocity.  A  slight 
decreases  in  the  slopes  of  both  (r>  and  (r,  >  are  also  ob¬ 
served.  From  Fig.  5.  the  mobilities  in  the  linear  region  are 
calculated  to  be'  8800.  520.  and  1 10  enr  ,  V  s  for  the  F.  /..  and 
.V  valleys,  respectively 

The  average  electron  energy  in  each  v  alley  relative  to  its 
lowest  energy  are  illustrated  in  Fig.  b.  For  low  fields,  the 
electron  energy  in  the  T  valley  increases  rapullv  with  field. 
In  tiiis  range  the  dominant  scattering  process  is  invariant 


polar  scattering.  Above  10  kV/cm.  where  population  inver¬ 
sion  is  observed  in  the  T  valley,  the  strong  intervalley  scat¬ 
tering  process  takes  energetic  electrons  out  of  the  T  valley, 
tints  reducing  the  rate  of  average  energy  increase.  In  the 
upper  v  alley  s,  the  average  energy  varies  almost  linearly  with 
field.  This  is  due  to  the  high  intervalley  scattering  rate  in  the 
upper  v  alleys. 

IV.  CONCLUDING  REMARKS 

In  this  paper,  a  technique  for  implementing  the  Monte 
Carlo  Method  with  self-scattering  has  been  presented.  This 
technique  leads  to  large  savings  in  computational  time,  and 
thus  allows,  for  a  given  CPU  time,  an  increase  in  the  number 
of  lest  particles  used  in  the  simulation.  This  reduces  the  sta¬ 
tistical  fluctuation  in  the  calculated  quantities.  With  a  code 
that  uses  this  technique,  we  have  investigated  the  steady- 
state  behavior  of  electrons  in  a  three-valley  model  of  GaAs 
and  influenced  by  a  dc  electric  field. 


FIG.  5.  (a)  Total  average  velocity  <t’>  and  the 
average  velocity  in  each  valley  as  a  function  of 
applied  electric  field.  Also  included  are  the  ex¬ 
perimental  results  by  Houston  and  Evans  (see 
Ref.  20),  Ruch  and  Kino  (see  Ref.  25),  and 
Riginos  (see  Ref.  26).  (b)  Total  average  ve¬ 
locity  ( v ),  the  average  velocity  in  T  valley 
(vr),  and  the  fraction  of  electrons  in  upper 
valleys  (short  dashed  line )  as  a  function  of  ap¬ 
plied  electric  field  for  fields  near  and  below  the 
velocity  maximum. 
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FIG.  6.  Mean  energy  of  electrons  in  each  valley 
relative  io  its  lowest  energy  as  a  function  of  ap¬ 
plied  electric  field. 


At  low  fields,  the  results  obtained  agree  well  with  pre- 
v  unis  theoretical  and  experimental  results.  At  high  fields,  we 
have  observed  a  second  transition  region  where  the  average 
velocity  in  the  F  valley  rapidly  changes  with  field.  This 
change  is  due  to  the  increase  in  the  population  of  the  X  val- 
ie\s  and  the  large  A-to-T  scattering  rate.  This  process  also 
has  an  effect  on  the  behavior  of  the  distribution  function, 
electron  population  in  each  valley,  and  the  transport  param¬ 
eters  for  fields  at  and  above  the  transition  region.  In  particu¬ 


lar.  the  population  inversion  observed  in  the  T  valley  for 
fields  above  10  kV/cm  is  not  as  pronounced  (as  it  is  in  a  two- 
valley  model). 

In  fast  transient  situations,  due  to  the  differences  in  scat¬ 
tering  rates,  a  three-valley  model  of  GaAs  is  desirable  for 
describing  the  dynamics  of  the  carriers.  The  results  present¬ 
ed  in  this  paper  can  be  used  in  the  development  of  a  three- 
valley  hydrodynamic  model. 
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a  macro-kinetic  equation1  which,  together  with  the  finite  set  of 
moment  equations,  constitutes  the  model.  This  set  can  be  used  to 
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times,  respectively.  The  first  (Si),  second  (S2)  and  third  (S3) 
levels  of  descriptions  are  valid  for  times  in  the  order  of  r,  re ,  and 
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are  presented. 
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